Show code
# Create two scenarios with different steepness values
scenario_low <- create_scenario(h = 0.6, "Low,h=0.6")
scenario_high <- create_scenario(h = 0.9, "High, h=0.9")
scenarios <- list(scenario_low, scenario_high)Comparing Low and High Steepness Scenarios with constant HCR
Recruitment productivity under climate change is uncertain, and stock-recruitment steepness (\(h\)) is a key parameter that drives sustainable harvest levels. Evidence from climate-ecosystem modeling and recruitment studies suggests steepness and productivity can shift with warming and prey dynamics (Hollowed et al. 2020; Holsman et al. 2020; Spencer et al. 2019; Szuwalski et al. 2023). Because steepness informs reference points and F proxies, mis-specification can change expected yields and risk profiles (Szuwalski and Punt 2025; Punt et al. 2024).
At the same time, there is pressure to incorporate environmental covariates into harvest control rules (HCRs) to improve adaptability. These approaches can help track productivity changes, but they also add interpretive uncertainty and require strong validation (Punt et al. 2024; Szuwalski et al. 2023). In contrast, tier-based HCRs remain transparent and easier to communicate, with clear links from biomass to fishing mortality (North Pacific Fishery Management Council 2024). This study asks whether a simple, transparent HCR can be robust to plausible productivity shifts without relying on covariates.
We use a stylized pollock-like, age-structured model and compare two steepness scenarios (h=0.6 and h=0.9) (North Pacific Fishery Management Council 2024). The HCR is held constant (F40% SPR with a biomass ramp from 0.4\(B_0\) to 0.05\(B_0\)), consistent with the Tier 3 framework used for Bering Sea groundfish (North Pacific Fishery Management Council 2024; Ianelli et al. 2011). The goal is to evaluate whether a single rule performs reasonably across both low and high productivity regimes, and to quantify the tradeoffs if productivity assumptions are wrong (Szuwalski and Punt 2025).
This two-scenario comparison also frames a plausible climate-driven trajectory: productivity could shift over time from a high-steepness regime to a low-steepness regime (or vice versa) as environmental conditions change. Interpreting the scenarios as endpoints of a time-varying process clarifies the management problem: the HCR is calibrated under one regime but must remain acceptable if the system transitions to the other. The results therefore inform both static mis-specification risk and the dynamic case of a climate-driven regime shift.
Key design choices include: - 15 age classes with Bering Sea-like growth and maturity schedules (Ianelli et al. 2011) - Beverton-Holt recruitment with lognormal variability (CV 0.7) (Spencer et al. 2019) - 50-year projections with 100 Monte Carlo replicates per scenario - No explicit environmental covariates; steepness shifts serve as a reduced-form proxy for climate-driven productivity changes (Hollowed et al. 2020; Holsman et al. 2020)
Outputs include reference points, time series of biomass and catch, and sensitivity comparisons across steepness values. We move from model setup to reference points in Table 1, then visualize recruitment and yield tradeoffs in Figure 1 and Figure 2. Simulation outcomes are summarized in Figure 4 and Figure 5, with sensitivity results in Figure 7, similar to HCR robustness evaluations in recent studies (Punt et al. 2024).
This section documents the equations and inputs used to generate the reference points and simulations shown later.
We define the annual accounting that drives the time series trajectories shown in Figure 4.
The population is tracked as numbers-at-age \(N_{a,t}\) where \(a\) denotes age and \(t\) denotes year. The population evolves according to:
Recruitment (age 1): \[N_{1,t} = R_t\]
Ages 2 through 14: \[N_{a,t} = N_{a-1,t-1} \cdot e^{-Z_{a-1,t-1}}\]
Plus group (age 15): \[N_{15,t} = N_{14,t-1} \cdot e^{-Z_{14,t-1}} + N_{15,t-1} \cdot e^{-Z_{15,t-1}}\]
where \(Z_{a,t} = M + F_t \cdot s_a\) is the total instantaneous mortality rate.
These equations determine the age structure that feeds into biomass and catch trajectories in Figure 4.
Biological inputs control growth, maturity, and vulnerability, which in turn shape the yield curves in Figure 2 and the SSB trends in Figure 4.
Weight-at-age is specified empirically based on the updated schedule:
| Age | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15+ |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Weight (kg) | 0.029 | 0.180 | 0.372 | 0.491 | 0.613 | 0.749 | 0.883 | 1.014 | 1.133 | 1.243 | 1.348 | 1.388 | 1.469 | 1.553 | 1.791 |
These weights translate numbers-at-age into biomass and catch, which are summarized in Figure 4 and Figure 2.
Maturity is specified directly using fixed proportions at age, while fishery selectivity follows a logistic function:
| Age | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15+ |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Proportion mature | 0.000 | 0.008 | 0.289 | 0.641 | 0.842 | 0.901 | 0.947 | 0.963 | 0.970 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 |
\[s_a = \frac{1}{1 + e^{-\gamma_s(a - a_{50}^s)}} \quad \text{with } a_{50}^s = 4.0\]
Maturity governs SSB while selectivity determines fishing impact, which together shape the trajectories in Figure 4.
SSB is defined as female spawning biomass only:
\[SSB_t = \sum_{a=1}^{15} N_{a,t} \cdot W_a \cdot m_a \cdot 0.5\]
This definition is used directly in the SSB time series shown in Figure 4 and the Monte Carlo envelope in Figure 5.
Recruitment dynamics are summarized visually in Figure 1 and influence the variability seen in Figure 5.
The expected recruitment follows the Beverton-Holt stock-recruitment relationship:
\[R_t = \frac{4hR_0 \cdot SSB_{t-1}}{SSB_0(1-h) + SSB_{t-1}(5h-1)}\]
Where steepness (\(h\)) represents the fraction of \(R_0\) produced when SSB is at 20% of unfished levels. Higher steepness indicates a more resilient stock that maintains recruitment even at lower biomass levels.
As shown in Figure 1, this formulation differs across the two steepness scenarios.
Recruitment includes multiplicative log-normal process error with CV = 0.7:
\[R_t = \bar{R}_t \cdot e^{\epsilon_t - \sigma_R^2/2}\]
where \(\epsilon_t \sim N(0, \sigma_R^2)\) and \(\sigma_R = \sqrt{\ln(1 + 0.7^2)} \approx 0.614\).
This stochasticity drives the spread in Figure 5 and contributes to the variability summarized in Table 4.
We instantiate the two steepness scenarios used throughout the tables and figures.
# Create two scenarios with different steepness values
scenario_low <- create_scenario(h = 0.6, "Low,h=0.6")
scenario_high <- create_scenario(h = 0.9, "High, h=0.9")
scenarios <- list(scenario_low, scenario_high)These objects feed the reference points in Table 1 and the graphical comparisons in Figure 1 and Figure 2.
All simulations initialize the population at the scenario’s Fmsy equilibrium rather than the unfished B0 state.
We summarize equilibrium reference points for each scenario to quantify differences in Fmsy, Bmsy, and MSY (see Table 1).
comparison_df <- compare_reference_points(scenarios)
knitr::kable(comparison_df, align = c("l", "c", "c", "c", "c", "c", "c", "c"))| Scenario | Steepness | Fmsy | Bmsy (Mt) | MSY (Mt) | SSB0 (Mt) | Bmsy/SSB0 | R_msy (B) |
|---|---|---|---|---|---|---|---|
| Low,h=0.6 | 0.6 | 0.3517 | 2.66 | 1.31 | 8.29 | 0.321 | 9.98 |
| High, h=0.9 | 0.9 | 1.2805 | 1.59 | 2.00 | 8.29 | 0.192 | 12.08 |
Table 1 provides the numeric basis for the contrasts discussed below and for the vertical markers in Figure 2.
ref_low <- scenario_low$ref_points
ref_high <- scenario_high$ref_pointsWe expand on the table summary by translating the reference point differences into percent changes.
The contrast between scenarios reveals important management implications:
| Metric | Low Steepness (h=0.6) | High Steepness (h=0.9) | Change |
|---|---|---|---|
| \(F_{MSY}\) | 0.352 | 1.28 | 264% higher |
| \(B_{MSY}\) (Mt) | 2.66 | 1.59 | -40% lower |
| MSY (Mt) | 1.31 | 2 | 53% higher |
| \(B_{MSY}/SSB_0\) | 0.321 | 0.192 | -40% lower |
Interpretation:
These figures visualize how the two steepness settings diverge in recruitment, yield, and harvest control (see Figure 1, Figure 2, and Figure 3).
As shown in Figure 1, the stock-recruitment curves differ and highlight the different Bmsy locations.
plot_sr_comparison(scenarios, colors = c("blue", "red"))The stock-recruitment curves illustrate how steepness affects the relationship between spawning biomass and recruitment:
As shown in Figure 2, equilibrium yield responds differently to fishing mortality across scenarios.
plot_yield_comparison(scenarios, colors = c("blue", "red"))The yield curves show:
We add a sloping harvest control rule based on \(F_{40\%}\) (SPR) and unfished biomass \(B_0\). The target F values are reported in Table 2 and the rule shape is shown in Figure 3.
\[ F_t = \begin{cases} F_{40\%} & \text{if } SSB_t \ge 0.4 B_0 \\ F_{40\%} \cdot \dfrac{SSB_t - 0.05 B_0}{0.4 B_0 - 0.05 B_0} & \text{if } 0.05 B_0 < SSB_t < 0.4 B_0 \\ 0 & \text{if } SSB_t \le 0.05 B_0 \end{cases} \]
F40_low <- calc_F_spr(scenario_low$params, spr_target = 0.4)
F40_high <- calc_F_spr(scenario_high$params, spr_target = 0.4)
f40_table <- data.frame(
Scenario = c("Low steepness (h=0.6)", "High steepness (h=0.9)"),
`F40%` = round(c(F40_low, F40_high), 3),
check.names = FALSE
)
knitr::kable(f40_table)| Scenario | F40% |
|---|---|
| Low steepness (h=0.6) | 0.412 |
| High steepness (h=0.9) | 0.412 |
ssb_ratio <- seq(0, 1, by = 0.01)
F_hcr_low <- calc_sloping_hcr(ssb_ratio * ref_low$SSB0, ref_low$SSB0, F40_low)
F_hcr_high <- calc_sloping_hcr(ssb_ratio * ref_high$SSB0, ref_high$SSB0, F40_high)
y_max <- max(c(F_hcr_low, F_hcr_high)) * 1.1
hcr_df <- data.frame(
SSB_ratio = rep(ssb_ratio, 2),
F = c(F_hcr_low, F_hcr_high),
Scenario = rep(c("h=0.6 (F40%)", "h=0.9 (F40%)"), each = length(ssb_ratio))
)
hcr_colors <- c("h=0.6 (F40%)" = "blue", "h=0.9 (F40%)" = "red")
p_hcr <- ggplot(hcr_df, aes(x = SSB_ratio, y = F, color = Scenario)) +
geom_line(size = 1) +
geom_vline(xintercept = c(0.05, 0.4), linetype = "dotted", color = "gray50") +
coord_cartesian(ylim = c(0, y_max)) +
labs(x = "SSB / B0", y = "Fishing Mortality (F)", title = "Harvest Control Rule") +
scale_color_manual(values = hcr_colors) +
ggthemes::theme_few()
p_hcrWe first compare single realizations at Fmsy (Figure 4), then summarize means (Table 3) and variability across Monte Carlo runs (Figure 5 and Table 4).
As shown in Figure 4, we plot a single 50-year trajectory for each scenario at its Fmsy.
# Run simulations
sim_low <- run_simulation(scenario_low$params, F_series = ref_low$Fmsy, seed = 42)
sim_high <- run_simulation(scenario_high$params, F_series = ref_high$Fmsy, seed = 42)
years <- 1:50
make_long <- function(sim, scenario_label) {
data.frame(
Year = rep(years, 4),
Scenario = scenario_label,
Metric = rep(c("SSB (Mt)", "Recruitment (B)", "Catch (Mt)", "F"), each = length(years)),
Value = c(sim$SSB / 1e9, sim$Recruitment / 1e9, sim$Catch / 1e9, sim$F)
)
}
sim_long <- rbind(
make_long(sim_low, "h=0.6"),
make_long(sim_high, "h=0.9")
)
sim_long$Metric <- factor(
sim_long$Metric,
levels = c("SSB (Mt)", "Recruitment (B)", "Catch (Mt)", "F")
)
ref_lines <- data.frame(
Metric = c("SSB (Mt)", "SSB (Mt)", "Catch (Mt)", "Catch (Mt)"),
Scenario = c("h=0.6", "h=0.9", "h=0.6", "h=0.9"),
Value = c(ref_low$Bmsy / 1e9, ref_high$Bmsy / 1e9, ref_low$MSY / 1e9, ref_high$MSY / 1e9)
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_single <- ggplot(sim_long, aes(x = Year, y = Value, color = Scenario)) +
geom_line(size = 1) +
geom_hline(
data = ref_lines,
aes(yintercept = Value, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE,
show.legend = FALSE
) +
facet_wrap(~ Metric, scales = "free_y", ncol = 2) +
scale_color_manual(values = scenario_colors) +
labs(x = "Year", y = NULL, title = "50-year simulations at Fmsy") +
ggthemes::theme_few()
p_singleTable 3 aggregates the single-simulation outcomes into mean and CV metrics.
summary_df <- data.frame(
Scenario = c("Low, h=0.6", "High, h=0.9"),
`Mean SSB (Mt)` = c(round(mean(sim_low$SSB)/1e9, 2), round(mean(sim_high$SSB)/1e9, 2)),
`CV SSB` = c(round(sd(sim_low$SSB)/mean(sim_low$SSB), 3),
round(sd(sim_high$SSB)/mean(sim_high$SSB), 3)),
`Mean Catch (Mt)` = c(round(mean(sim_low$Catch)/1e9, 2), round(mean(sim_high$Catch)/1e9, 2)),
`Mean Rec (B)` = c(round(mean(sim_low$Recruitment)/1e9, 2),
round(mean(sim_high$Recruitment)/1e9, 2)),
check.names = FALSE
)
knitr::kable(summary_df)| Scenario | Mean SSB (Mt) | CV SSB | Mean Catch (Mt) | Mean Rec (B) |
|---|---|---|---|---|
| Low, h=0.6 | 2.75 | 0.352 | 1.35 | 10.09 |
| High, h=0.9 | 1.60 | 0.413 | 2.01 | 12.13 |
As shown in Figure 5, stochastic variability across 100 runs is summarized with metrics reported in Table 4.
# Run Monte Carlo simulations
n_sims <- 100
n_years <- 50
ssb_low <- matrix(NA, n_years, n_sims)
ssb_high <- matrix(NA, n_years, n_sims)
catch_low <- matrix(NA, n_years, n_sims)
catch_high <- matrix(NA, n_years, n_sims)
for (i in 1:n_sims) {
sim_l <- run_simulation(scenario_low$params, F_series = ref_low$Fmsy, seed = i)
sim_h <- run_simulation(scenario_high$params, F_series = ref_high$Fmsy, seed = i)
ssb_low[, i] <- sim_l$SSB
ssb_high[, i] <- sim_h$SSB
catch_low[, i] <- sim_l$Catch
catch_high[, i] <- sim_h$Catch
}
ssb_low_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
SSB = as.vector(ssb_low) / 1e9,
Scenario = "h=0.6"
)
ssb_high_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
SSB = as.vector(ssb_high) / 1e9,
Scenario = "h=0.9"
)
ssb_df <- rbind(ssb_low_df, ssb_high_df)
median_df <- rbind(
data.frame(Year = 1:n_years, SSB = apply(ssb_low, 1, median) / 1e9, Scenario = "h=0.6"),
data.frame(Year = 1:n_years, SSB = apply(ssb_high, 1, median) / 1e9, Scenario = "h=0.9")
)
bmsy_df <- data.frame(
Scenario = c("h=0.6", "h=0.9"),
Bmsy = c(ref_low$Bmsy, ref_high$Bmsy) / 1e9
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_mc <- ggplot(ssb_df, aes(x = Year, y = SSB, group = interaction(Scenario, Sim), color = Scenario)) +
geom_line(alpha = 0.1, size = 0.4, show.legend = FALSE) +
geom_line(
data = median_df,
aes(x = Year, y = SSB, color = Scenario),
size = 1.1,
inherit.aes = FALSE
) +
geom_hline(
data = bmsy_df,
aes(yintercept = Bmsy, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE
) +
facet_wrap(~ Scenario, ncol = 2, scales = "free_y") +
scale_color_manual(values = scenario_colors) +
labs(x = "Year", y = "Female SSB (Mt)", title = "Monte Carlo SSB trajectories") +
ggthemes::theme_few() +
guides(color = "none")
p_mccatch_low_hcr <- matrix(NA, n_years, n_sims)
catch_high_hcr <- matrix(NA, n_years, n_sims)
for (i in 1:n_sims) {
sim_l_hcr <- run_simulation_hcr(
scenario_low$params,
F_target = F40_low,
B0 = ref_low$SSB0,
catch_cap = NULL,
seed = 2000 + i
)
sim_h_hcr <- run_simulation_hcr(
scenario_high$params,
F_target = F40_high,
B0 = ref_high$SSB0,
catch_cap = NULL,
seed = 3000 + i
)
catch_low_hcr[, i] <- sim_l_hcr$Catch
catch_high_hcr[, i] <- sim_h_hcr$Catch
}
catch_low_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_low) / 1e9,
Scenario = "h=0.6",
Policy = "Const. F"
)
catch_high_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_high) / 1e9,
Scenario = "h=0.9",
Policy = "Const. F"
)
catch_low_hcr_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_low_hcr) / 1e9,
Scenario = "h=0.6",
Policy = "FSPR HCR"
)
catch_high_hcr_df <- data.frame(
Year = rep(1:n_years, times = n_sims),
Sim = rep(1:n_sims, each = n_years),
Catch = as.vector(catch_high_hcr) / 1e9,
Scenario = "h=0.9",
Policy = "FSPR HCR"
)
catch_df <- rbind(catch_low_df, catch_high_df, catch_low_hcr_df, catch_high_hcr_df)
median_catch_df <- rbind(
data.frame(
Year = 1:n_years,
Catch = apply(catch_low, 1, median) / 1e9,
Scenario = "h=0.6",
Policy = "Const. F"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_high, 1, median) / 1e9,
Scenario = "h=0.9",
Policy = "Const. F"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_low_hcr, 1, median) / 1e9,
Scenario = "h=0.6",
Policy = "FSPR HCR"
),
data.frame(
Year = 1:n_years,
Catch = apply(catch_high_hcr, 1, median) / 1e9,
Scenario = "h=0.9",
Policy = "FSPR HCR"
)
)
policy_colors <- c("Const. F" = "#1F77B4", "FSPR HCR" = "#D95F02")
p_catch <- ggplot(
catch_df,
aes(x = Year, y = Catch, group = interaction(Scenario, Policy, Sim), color = Policy)
) +
geom_line(alpha = 0.08, size = 0.4, show.legend = FALSE) +
geom_line(
data = median_catch_df,
aes(x = Year, y = Catch, color = Policy),
size = 1.1,
inherit.aes = FALSE
) +
facet_wrap(~ Scenario, ncol = 2) +
scale_color_manual(values = policy_colors) +
labs(x = "Year", y = "Catch (Mt)", title = "Monte Carlo catch trajectories") +
ggthemes::theme_few()
p_catchmc_stats <- data.frame(
Scenario = c("Low, h=0.6", "High, h=0.9"),
`Mean terminal SSB (Mt)` = c(round(mean(ssb_low[n_years,])/1e9, 2),
round(mean(ssb_high[n_years,])/1e9, 2)),
`P(SSB < Bmsy) terminal` = c(round(mean(ssb_low[n_years,] < ref_low$Bmsy), 2),
round(mean(ssb_high[n_years,] < ref_high$Bmsy), 2)),
`Mean catch (Mt)` = c(round(mean(catch_low)/1e9, 2),
round(mean(catch_high)/1e9, 2)),
`CV mean catch` = c(round(sd(colMeans(catch_low))/mean(catch_low), 3),
round(sd(colMeans(catch_high))/mean(catch_high), 3)),
check.names = FALSE
)
knitr::kable(mc_stats)| Scenario | Mean terminal SSB (Mt) | P(SSB < Bmsy) terminal | Mean catch (Mt) | CV mean catch |
|---|---|---|---|---|
| Low, h=0.6 | 2.58 | 0.56 | 1.26 | 0.134 |
| High, h=0.9 | 1.64 | 0.54 | 1.96 | 0.105 |
To stress-test harvest rates, we run each scenario at constant fishing mortality values spanning 0.5-1.5 x \(F_{MSY}\). Table 5 lists the constant-F results alongside the FSPR HCR outcomes (with and without the 1.45 Mt cap), and Figure 7 visualizes the tradeoffs.
F_multipliers <- c(0.5, 0.75, 1.0, 1.25, 1.5)
n_sims_F <- 100
f_low <- run_F_sensitivity(scenario_low, F_multipliers, n_sims = n_sims_F, seed = 200)
f_high <- run_F_sensitivity(scenario_high, F_multipliers, n_sims = n_sims_F, seed = 500)
f_summary <- rbind(f_low, f_high)
f_summary$Policy <- "Const. F"
run_hcr_summary_table <- function(scenario, F_target, B0, n_sims, seed_base, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
ssb_mean <- numeric(n_sims)
ssb_terminal <- numeric(n_sims)
catch_mean <- numeric(n_sims)
f_mean <- numeric(n_sims)
ssb_lt_bmsy <- logical(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
ssb_mean[i] <- mean(sim$SSB)
ssb_terminal[i] <- sim$SSB[n_years]
catch_mean[i] <- mean(sim$Catch)
f_mean[i] <- mean(sim$F)
ssb_lt_bmsy[i] <- sim$SSB[n_years] < ref_points$Bmsy
}
data.frame(
Scenario = scenario$name,
Steepness = params$h,
Policy = policy_label,
F_multiplier = mean(f_mean) / ref_points$Fmsy,
F_applied = mean(f_mean),
Mean_SSB = mean(ssb_mean),
Terminal_SSB = mean(ssb_terminal),
P_SSB_lt_Bmsy = mean(ssb_lt_bmsy),
Mean_Catch = mean(catch_mean)
)
}
catch_cap <- 1.45e9
hcr_low <- run_hcr_summary_table(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 800,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high <- run_hcr_summary_table(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 900,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap <- run_hcr_summary_table(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 1000,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap <- run_hcr_summary_table(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 1100,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_summary <- rbind(hcr_low, hcr_high, hcr_low_cap, hcr_high_cap)
f_summary_all <- rbind(f_summary, hcr_summary)
f_table <- data.frame(
Scenario = f_summary_all$Scenario,
Policy = f_summary_all$Policy,
`F multiplier` = signif(f_summary_all$F_multiplier, 2),
`F applied` = round(f_summary_all$F_applied, 3),
`Mean SSB (Mt)` = round(f_summary_all$Mean_SSB / 1e9, 2),
`Terminal SSB (Mt)` = round(f_summary_all$Terminal_SSB / 1e9, 2),
`P(SSB < Bmsy) terminal` = round(f_summary_all$P_SSB_lt_Bmsy, 2),
`Mean Catch (Mt)` = round(f_summary_all$Mean_Catch / 1e9, 2),
check.names = FALSE
)
knitr::kable(f_table)| Scenario | Policy | F multiplier | F applied | Mean SSB (Mt) | Terminal SSB (Mt) | P(SSB < Bmsy) terminal | Mean Catch (Mt) |
|---|---|---|---|---|---|---|---|
| Low,h=0.6 | Const. F | 0.50 | 0.176 | 4.11 | 4.32 | 0.08 | 1.10 |
| Low,h=0.6 | Const. F | 0.75 | 0.264 | 3.22 | 3.30 | 0.29 | 1.24 |
| Low,h=0.6 | Const. F | 1.00 | 0.352 | 2.67 | 2.50 | 0.61 | 1.31 |
| Low,h=0.6 | Const. F | 1.20 | 0.440 | 2.19 | 2.17 | 0.76 | 1.29 |
| Low,h=0.6 | Const. F | 1.50 | 0.528 | 1.87 | 1.72 | 0.90 | 1.27 |
| High, h=0.9 | Const. F | 0.50 | 0.640 | 2.37 | 2.52 | 0.08 | 1.85 |
| High, h=0.9 | Const. F | 0.75 | 0.960 | 1.85 | 1.90 | 0.45 | 1.93 |
| High, h=0.9 | Const. F | 1.00 | 1.280 | 1.55 | 1.45 | 0.67 | 1.95 |
| High, h=0.9 | Const. F | 1.20 | 1.601 | 1.39 | 1.29 | 0.77 | 2.01 |
| High, h=0.9 | Const. F | 1.50 | 1.921 | 1.21 | 1.19 | 0.89 | 1.96 |
| Low,h=0.6 | FSPR HCR | 0.90 | 0.318 | 2.72 | 2.69 | 0.57 | 1.26 |
| High, h=0.9 | FSPR HCR | 0.28 | 0.354 | 3.22 | 3.25 | 0.00 | 1.63 |
| Low,h=0.6 | FSPR HCR cap | 0.78 | 0.275 | 3.13 | 3.20 | 0.35 | 1.21 |
| High, h=0.9 | FSPR HCR cap | 0.19 | 0.237 | 4.22 | 4.40 | 0.00 | 1.35 |
plot_data <- f_summary
plot_data$Policy <- "Const. F"
plot_data$Scenario_plot <- ifelse(
plot_data$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
catch_cap <- 1.45e9
run_hcr_summary <- function(scenario, F_target, B0, n_sims, seed_base, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
ssb_mean <- numeric(n_sims)
ssb_terminal <- numeric(n_sims)
catch_mean <- numeric(n_sims)
f_mean <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
ssb_mean[i] <- mean(sim$SSB)
ssb_terminal[i] <- sim$SSB[n_years]
catch_mean[i] <- mean(sim$Catch)
f_mean[i] <- mean(sim$F)
}
data.frame(
Scenario = scenario$name,
F_multiplier = mean(f_mean) / ref_points$Fmsy,
Mean_SSB = mean(ssb_mean),
Terminal_SSB = mean(ssb_terminal),
Mean_Catch = mean(catch_mean),
Policy = policy_label
)
}
hcr_low <- run_hcr_summary(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 800,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high <- run_hcr_summary(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 900,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap <- run_hcr_summary(
scenario_low,
F40_low,
ref_low$SSB0,
n_sims_F,
seed_base = 1000,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap <- run_hcr_summary(
scenario_high,
F40_high,
ref_high$SSB0,
n_sims_F,
seed_base = 1100,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_summary <- rbind(hcr_low, hcr_high, hcr_low_cap, hcr_high_cap)
hcr_summary$Scenario_plot <- ifelse(
hcr_summary$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
plot_long <- rbind(
data.frame(
F_multiplier = plot_data$F_multiplier,
Scenario = plot_data$Scenario_plot,
Metric = "Mean Catch (Mt)",
Value = plot_data$Mean_Catch / 1e9,
Policy = plot_data$Policy
),
data.frame(
F_multiplier = plot_data$F_multiplier,
Scenario = plot_data$Scenario_plot,
Metric = "Terminal SSB (Mt)",
Value = plot_data$Terminal_SSB / 1e9,
Policy = plot_data$Policy
),
data.frame(
F_multiplier = hcr_summary$F_multiplier,
Scenario = hcr_summary$Scenario_plot,
Metric = "Mean Catch (Mt)",
Value = hcr_summary$Mean_Catch / 1e9,
Policy = hcr_summary$Policy
),
data.frame(
F_multiplier = hcr_summary$F_multiplier,
Scenario = hcr_summary$Scenario_plot,
Metric = "Terminal SSB (Mt)",
Value = hcr_summary$Terminal_SSB / 1e9,
Policy = hcr_summary$Policy
)
)
plot_long$Metric <- factor(plot_long$Metric, levels = c("Mean Catch (Mt)", "Terminal SSB (Mt)"))
bmsy_df <- data.frame(
Scenario = c("h=0.6", "h=0.9"),
Metric = "Terminal SSB (Mt)",
Value = c(ref_low$Bmsy, ref_high$Bmsy) / 1e9
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_f <- ggplot() +
geom_line(
data = plot_long[plot_long$Policy == "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario),
size = 1
) +
geom_point(
data = plot_long[plot_long$Policy == "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario),
size = 2
) +
geom_point(
data = plot_long[plot_long$Policy != "Const. F", ],
aes(x = F_multiplier, y = Value, color = Scenario, shape = Policy),
size = 3,
stroke = 0.8
) +
geom_hline(
data = bmsy_df,
aes(yintercept = Value, color = Scenario),
linetype = "dashed",
size = 0.8,
inherit.aes = FALSE,
show.legend = FALSE
) +
facet_wrap(~ Metric, scales = "free_y", ncol = 2) +
scale_color_manual(values = scenario_colors) +
scale_shape_manual(values = c("FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "F multiplier (relative to Fmsy)", y = NULL, title = "F sensitivity tradeoffs") +
ggthemes::theme_few()
p_fn_sims_points <- n_sims_F
catch_cap <- 1.45e9
sim_constant_points <- function(scenario, f_multiplier, seed_base, n_sims) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
mean_catch <- numeric(n_sims)
terminal_ssb <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation(params, F_series = ref_points$Fmsy * f_multiplier, seed = seed_base + i)
mean_catch[i] <- mean(sim$Catch)
terminal_ssb[i] <- sim$SSB[n_years]
}
data.frame(
Scenario = scenario$name,
F_multiplier = rep(f_multiplier, n_sims),
Mean_Catch = mean_catch,
Terminal_SSB = terminal_ssb,
Policy = "Const. F"
)
}
sim_hcr_points <- function(scenario, F_target, B0, seed_base, n_sims, catch_cap, policy_label) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
mean_catch <- numeric(n_sims)
terminal_ssb <- numeric(n_sims)
mean_f <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
mean_catch[i] <- mean(sim$Catch)
terminal_ssb[i] <- sim$SSB[n_years]
mean_f[i] <- mean(sim$F)
}
data.frame(
Scenario = scenario$name,
F_multiplier = mean_f / ref_points$Fmsy,
Mean_Catch = mean_catch,
Terminal_SSB = terminal_ssb,
Policy = policy_label
)
}
const_low_pts <- sim_constant_points(scenario_low, 1.0, seed_base = 1200, n_sims_points)
const_high_pts <- sim_constant_points(scenario_high, 1.0, seed_base = 1300, n_sims_points)
hcr_low_pts <- sim_hcr_points(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1400,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high_pts <- sim_hcr_points(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1500,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap_pts <- sim_hcr_points(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1600,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap_pts <- sim_hcr_points(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1700,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
points_df <- rbind(const_low_pts, const_high_pts, hcr_low_pts, hcr_high_pts, hcr_low_cap_pts, hcr_high_cap_pts)
points_df$Scenario_plot <- ifelse(
points_df$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
bmsy_lookup <- ifelse(
points_df$Scenario_plot == "h=0.6",
ref_low$Bmsy,
ref_high$Bmsy
)
plot_points <- data.frame(
Catch_Mt = points_df$Mean_Catch / 1e9,
Terminal_SSB_Bmsy = points_df$Terminal_SSB / bmsy_lookup,
Scenario = points_df$Scenario_plot,
Policy = points_df$Policy,
PolicyGroup = ifelse(points_df$Policy == "Const. F", "Const. F", "SPR-based HCR")
)
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
p_fp <- ggplot(plot_points, aes(x = Catch_Mt, y = Terminal_SSB_Bmsy, color = Scenario, shape = Policy)) +
geom_point(alpha = 0.6, size = 2) +
geom_smooth(
data = plot_points,
aes(x = Catch_Mt, y = Terminal_SSB_Bmsy, linetype = PolicyGroup, group = PolicyGroup),
method = "loess",
se = FALSE,
size = 0.8,
color = "black",
inherit.aes = FALSE
) +
expand_limits(y = 0) +
scale_color_manual(values = scenario_colors) +
scale_linetype_manual(values = c("Const. F" = "solid", "SPR-based HCR" = "dashed")) +
scale_shape_manual(values = c("Const. F" = 16, "FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "Mean Catch (Mt)", y = "Terminal SSB / Bmsy", title = "F sensitivity: individual simulation draws") +
ggthemes::theme_few()
p_fpcatch_cap <- 1.45e9
last10_stats_constant <- function(scenario, f_multiplier, seed_base, n_sims) {
params <- scenario$params
ref_points <- scenario$ref_points
n_years <- params$n_years
last_years <- (n_years - 9):n_years
mean_catch <- numeric(n_sims)
cv_catch <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation(params, F_series = ref_points$Fmsy * f_multiplier, seed = seed_base + i)
catch_last <- sim$Catch[last_years]
mean_catch[i] <- mean(catch_last)
cv_catch[i] <- sd(catch_last) / mean(catch_last)
}
data.frame(
Scenario = scenario$name,
Mean_Catch = mean_catch,
Catch_CV = cv_catch,
Policy = "Const. F"
)
}
last10_stats_hcr <- function(scenario, F_target, B0, seed_base, n_sims, catch_cap, policy_label) {
params <- scenario$params
n_years <- params$n_years
last_years <- (n_years - 9):n_years
mean_catch <- numeric(n_sims)
cv_catch <- numeric(n_sims)
for (i in 1:n_sims) {
sim <- run_simulation_hcr(
params,
F_target = F_target,
B0 = B0,
catch_cap = catch_cap,
seed = seed_base + i
)
catch_last <- sim$Catch[last_years]
mean_catch[i] <- mean(catch_last)
cv_catch[i] <- sd(catch_last) / mean(catch_last)
}
data.frame(
Scenario = scenario$name,
Mean_Catch = mean_catch,
Catch_CV = cv_catch,
Policy = policy_label
)
}
const_low_var <- last10_stats_constant(scenario_low, 1.0, seed_base = 1600, n_sims_points)
const_high_var <- last10_stats_constant(scenario_high, 1.0, seed_base = 1700, n_sims_points)
hcr_low_var <- last10_stats_hcr(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 1800,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_high_var <- last10_stats_hcr(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 1900,
n_sims_points,
catch_cap = NULL,
policy_label = "FSPR HCR"
)
hcr_low_cap_var <- last10_stats_hcr(
scenario_low,
F40_low,
ref_low$SSB0,
seed_base = 2000,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
hcr_high_cap_var <- last10_stats_hcr(
scenario_high,
F40_high,
ref_high$SSB0,
seed_base = 2100,
n_sims_points,
catch_cap = catch_cap,
policy_label = "FSPR HCR cap"
)
var_df <- rbind(const_low_var, const_high_var, hcr_low_var, hcr_high_var, hcr_low_cap_var, hcr_high_cap_var)
var_df$Scenario_plot <- ifelse(
var_df$Scenario == scenario_low$name,
"h=0.6",
"h=0.9"
)
var_df$Mean_Catch_Mt <- var_df$Mean_Catch / 1e9
scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")
var_summary <- aggregate(
cbind(Mean_Catch_Mt, Catch_CV) ~ Scenario_plot + Policy,
data = var_df,
FUN = mean
)
p_var <- ggplot(var_summary, aes(x = Mean_Catch_Mt, y = Catch_CV, color = Scenario_plot, shape = Policy)) +
geom_point(size = 5) +
expand_limits(y = 0) +
scale_color_manual(values = scenario_colors) +
scale_shape_manual(values = c("Const. F" = 16, "FSPR HCR" = 17, "FSPR HCR cap" = 15)) +
labs(x = "Mean Catch (Mt, last 10 years)", y = "Catch CV (last 10 years)", title = "Catch variability vs mean catch") +
ggthemes::theme_few()
p_varTo represent a climate-driven regime shift, we run a separate simulation where steepness alternates between high and low values every 10 years. The capped SPR HCR (F40% with a 1.45 Mt catch cap) is applied throughout, and we summarize productivity using log(recruitment/SSB) for each simulation run.
n_years_alt <- scenario_low$params$n_years
block_years <- 10
block_id <- rep(seq_len(ceiling(n_years_alt / block_years)), each = block_years)[1:n_years_alt]
h_series <- ifelse(block_id %% 2 == 1, scenario_high$params$h, scenario_low$params$h)
n_sims_alt <- 100
catch_cap_alt <- 1.45e9
F_target_alt <- calc_F_spr(scenario_low$params, spr_target = 0.4)
B0_alt <- calc_SSB0(scenario_low$params)
log_rps_list <- vector("list", n_sims_alt)
for (i in 1:n_sims_alt) {
sim_alt <- run_simulation_hcr_variable_h(
scenario_low$params,
h_series = h_series,
F_target = F_target_alt,
B0 = B0_alt,
catch_cap = catch_cap_alt,
seed = 4000 + i
)
rps <- sim_alt$Recruitment[-1] / pmax(sim_alt$SSB[-length(sim_alt$SSB)], 1e-12)
log_rps_list[[i]] <- data.frame(
Year = 2:n_years_alt,
Sim = i,
Log_R_SSB = log(rps)
)
}
log_rps_df <- do.call(rbind, log_rps_list)
median_log_rps <- aggregate(Log_R_SSB ~ Year, data = log_rps_df, FUN = median)
p_alt <- ggplot(log_rps_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(x = "Year", y = "log(Recruitment / SSB)", title = "Alternating steepness (10-year blocks)") +
ggthemes::theme_few()
p_altWe repeat the alternating-steepness experiment with low recruitment variability (\(\sigma_R = 0.1\)) to show how reduced process error changes the log(recruitment/SSB) dynamics under the capped HCR.
n_years_alt_low <- scenario_low$params$n_years
block_years_low <- 10
block_id_low <- rep(seq_len(ceiling(n_years_alt_low / block_years_low)), each = block_years_low)[1:n_years_alt_low]
h_series_low <- ifelse(block_id_low %% 2 == 1, scenario_high$params$h, scenario_low$params$h)
params_low_sigma <- scenario_low$params
params_low_sigma$rec_sigma <- 0.1
params_low_sigma$rec_cv <- sqrt(exp(params_low_sigma$rec_sigma^2) - 1)
n_sims_alt_low <- 100
catch_cap_alt_low <- 1.45e9
F_target_alt_low <- calc_F_spr(params_low_sigma, spr_target = 0.4)
B0_alt_low <- calc_SSB0(params_low_sigma)
log_rps_list_low <- vector("list", n_sims_alt_low)
for (i in 1:n_sims_alt_low) {
sim_alt_low <- run_simulation_hcr_variable_h(
params_low_sigma,
h_series = h_series_low,
F_target = F_target_alt_low,
B0 = B0_alt_low,
catch_cap = catch_cap_alt_low,
seed = 5000 + i
)
rps_low <- sim_alt_low$Recruitment[-1] / pmax(sim_alt_low$SSB[-length(sim_alt_low$SSB)], 1e-12)
log_rps_list_low[[i]] <- data.frame(
Year = 2:n_years_alt_low,
Sim = i,
Log_R_SSB = log(rps_low)
)
}
log_rps_low_df <- do.call(rbind, log_rps_list_low)
median_log_rps_low <- aggregate(Log_R_SSB ~ Year, data = log_rps_low_df, FUN = median)
p_alt_low <- ggplot(log_rps_low_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps_low,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(
x = "Year",
y = "log(Recruitment / SSB)",
title = "Alternating steepness (10-year blocks), sigma_R = 0.1"
) +
ggthemes::theme_few()
p_alt_lowAs a contrast to the capped HCR, we apply the “correct” \(F_{MSY}\) for each steepness block (high or low) while keeping \(\sigma_R = 0.1\).
params_low_sigma_low <- params_low_sigma
params_low_sigma_low$h <- scenario_low$params$h
params_low_sigma_high <- params_low_sigma
params_low_sigma_high$h <- scenario_high$params$h
Fmsy_low_sigma <- calc_reference_points(params_low_sigma_low)$Fmsy
Fmsy_high_sigma <- calc_reference_points(params_low_sigma_high)$Fmsy
F_series_alt <- ifelse(block_id_low %% 2 == 1, Fmsy_high_sigma, Fmsy_low_sigma)
log_rps_list_fmsy <- vector("list", n_sims_alt_low)
for (i in 1:n_sims_alt_low) {
sim_alt_fmsy <- run_simulation_variable_h(
params_low_sigma,
h_series = h_series_low,
F_series = F_series_alt,
seed = 6000 + i
)
rps_fmsy <- sim_alt_fmsy$Recruitment[-1] / pmax(sim_alt_fmsy$SSB[-length(sim_alt_fmsy$SSB)], 1e-12)
log_rps_list_fmsy[[i]] <- data.frame(
Year = 2:n_years_alt_low,
Sim = i,
Log_R_SSB = log(rps_fmsy)
)
}
log_rps_fmsy_df <- do.call(rbind, log_rps_list_fmsy)
median_log_rps_fmsy <- aggregate(Log_R_SSB ~ Year, data = log_rps_fmsy_df, FUN = median)
p_alt_fmsy <- ggplot(log_rps_fmsy_df, aes(x = Year, y = Log_R_SSB, group = Sim)) +
geom_line(alpha = 0.2, size = 0.4, color = "steelblue") +
geom_line(
data = median_log_rps_fmsy,
aes(x = Year, y = Log_R_SSB),
inherit.aes = FALSE,
color = "black",
size = 1.1
) +
labs(
x = "Year",
y = "log(Recruitment / SSB)",
title = "Alternating steepness, regime-specific Fmsy (sigma_R = 0.1)"
) +
ggthemes::theme_few()
p_alt_fmsyEven with a known alternation in steepness, detecting a productivity shift in practice can be difficult. Recruitment variability, observation error in SSB, and the lagged relationship between SSB and recruitment can mask regime changes, especially over short time windows. This makes it challenging to distinguish a true shift in productivity from stochastic swings, and underscores the value of trigger rules that rely on sustained signals rather than single-year anomalies.
This comparison of low (\(h = 0.6\)) and high (\(h = 0.9\)) steepness scenarios demonstrates the critical influence of stock-recruitment steepness on fishery management. The main evidence comes from the reference points in Table 1 and the tradeoff patterns in Figure 1, Figure 2, and Figure 4. This simulation is a simplified, pollock-like evaluation of SPR-based control rules, grounded in the Tier 3 framework described in the BSAI FMP (North Pacific Fishery Management Council 2024). In this model, the proxy \(F_{MSY}\) is \(F_{40\%}\), with a sloping rule that holds \(F_{40\%}\) above 0.4 \(B_0\) and linearly ramps to zero at 0.05 \(B_0\). That structure aligns with SPR proxy guidance and the literature on SPR-based reference points (North Pacific Fishery Management Council 2024; Szuwalski and Punt 2025), while remaining intentionally stylized to match the scenarios tested here.
The two steepness scenarios (\(h = 0.6\) and \(h = 0.9\)) serve as a tractable proxy for productivity shifts discussed in climate and recruitment literature (Spencer et al. 2019; Punt et al. 2024; Szuwalski et al. 2023). We do not explicitly model environmental covariates or the ACLIM framework, but those references motivate the sensitivity analysis and the interpretation of performance differences across productivity regimes (Hollowed et al. 2020).
Finally, the catch-capped HCR variant (1.45 Mt) is a stylized analogue to ecosystem-cap effects that can dampen catch variability and conserve biomass during high abundance (Holsman et al. 2020; Ianelli et al. 2011). While the cap value here is illustrative and not a direct re-creation of the 2 Mt cap, the direction and intent are consistent with the cited findings, making the comparison appropriate for the simplified model context.
This analysis can be extended to explicitly test management objectives and triggers. Objectives define the desired outcomes and can be evaluated directly with the performance metrics already computed here (mean catch, catch variability, and risk of low biomass). Triggers then specify when the baseline capped HCR should tighten or relax specific elements of the rule (e.g., a lower F target, a smaller cap, or a steeper ramp).
Examples of objectives that map cleanly onto the current outputs include:
Illustrative triggers that could be implemented within the capped HCR framework include:
These triggers can be evaluated by their false-alarm rate, time-to-recovery, and tradeoffs against yield, allowing the Council to align management actions with explicit objectives rather than relying on implicit assumptions about productivity.
Illustrative trigger table (capped HCR adjustments):
| Trigger | Metric | Threshold | Capped HCR adjustment |
|---|---|---|---|
| Biomass safeguard | 3-year mean \(SSB/B_{MSY}\) | < 1.0 | Reduce \(F_{40\%}\) target by 20% for 5 years |
| Persistent decline | SSB trend | Declines in 4 of 5 years | Lower cap by 15% and steepen the ramp |
| Catch instability | Catch CV (last 10 years) | > 0.30 | Temporarily tighten cap by 10% |
| Recruitment anomaly | Recruitment percentile | < 20th percentile for 3 of 5 years | Reduce \(F_{40\%}\) by 15% and tighten ramp |
| Prolonged low recruitment | Recruitment vs mean | Below mean for 5 consecutive years | Reduce \(F_{40\%}\) by 10% and lower cap by 10% |
| Metric | h = 0.6 | h = 0.9 |
|---|---|---|
| \(F_{MSY}\) | 0.352 | 1.28 |
| \(B_{MSY}\) | 2.66 Mt | 1.59 Mt |
| MSY | 1.31 Mt | 2 Mt |
| \(B_{MSY}/SSB_0\) | 0.32 | 0.19 |
Higher steepness dramatically increases sustainable fishing mortality: \(F_{MSY}\) increases from 0.35 to 1.28 (a 264% increase).
MSY is 53% higher with high steepness (2 vs 1.31 million tonnes).
\(B_{MSY}\) is lower with high steepness because recruitment remains high even at reduced biomass levels.
Management implications: Uncertainty in steepness translates directly to uncertainty in sustainable harvest levels. A stock managed as if \(h = 0.9\) when true \(h = 0.6\) would be at high risk of overfishing.
F40% sloping HCR vs constant F: The proxy Fmsy (F40%) with the sloping control rule yields a different catch/SSB tradeoff than the constant-F multipliers, and its mean realized F/Fmsy and outcomes shift between \(h = 0.6\) and \(h = 0.9\) (see Figure 7 and Figure 8).
Session information is reported below to document the computational environment used to generate the figures and tables.
sessionInfo()R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggthemes_5.2.0 ggplot2_4.0.1
loaded via a namespace (and not attached):
[1] Matrix_1.7-4 gtable_0.3.6 jsonlite_2.0.0 dplyr_1.1.4
[5] compiler_4.5.2 tidyselect_1.2.1 stringr_1.6.0 splines_4.5.2
[9] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 lattice_0.22-7
[13] R6_2.6.1 labeling_0.4.3 generics_0.1.4 knitr_1.50
[17] htmlwidgets_1.6.4 tibble_3.3.0 pillar_1.11.1 RColorBrewer_1.1-3
[21] rlang_1.1.6 stringi_1.8.7 xfun_0.55 S7_0.2.1
[25] cli_3.6.5 withr_3.0.2 magrittr_2.0.4 mgcv_1.9-3
[29] digest_0.6.39 grid_4.5.2 nlme_3.1-168 lifecycle_1.0.4
[33] vctrs_0.6.5 evaluate_1.0.5 glue_1.8.0 farver_2.1.2
[37] codetools_0.2-20 rmarkdown_2.30 purrr_1.2.0 tools_4.5.2
[41] pkgconfig_2.0.3 htmltools_0.5.9